# libraries
library(devtools)
install_github("htx-r/doseresNMA",force=TRUE)
library(doseresNMA)
library(R2jags)
library(rms)
library(meta)
library('MBNMAdose')
#** Prepare data
# load
mydata <- read.csv('DOSEmainanalysis.csv')
antidepdata0=mydata[mydata$exc==F,]
antidepdata <- antidepdata0[antidepdata0$Study_No!=161001,]
length(table(antidepdata0$Dose_delivered_mean))
# add OR to the dataset
antidepdata$studyid <- as.numeric(as.factor(antidepdata$Study_No))
antidepdata$nonResponders <- antidepdata$No_randomised- antidepdata$Responders
logORmat <- sapply(unique(antidepdata$studyid),function(i) createORreference.fun(antidepdata$Responders[antidepdata$studyid==i],antidepdata$No_randomised[antidepdata$studyid==i]),simplify = FALSE)
logORmat <- do.call(rbind,logORmat)
antidepdata$logOR <- c(logORmat[,1])
antidepdata$selogOR <- c(logORmat[,2])
#** JAGS data format
class <- as.numeric(as.factor(antidepdata$Drug))
class[class==5] <- 'A'
class[class==1|class==2|class==3] <- 'B'
class[class==4|class==6] <- 'C'
class <- as.numeric(as.factor(class))
mydata <- data.frame(studyid=antidepdata$Study_No,drug=antidepdata$Drug,dose=antidepdata$Dose_delivered_mean,
r=antidepdata$Responders,n=antidepdata$No_randomised)
# ROB=antidepdata$Overall_study_RoB,Study_year=antidepdata$Study_year,selogOR=antidepdata$selogOR,class=class)
dosresnetmeta.jagsdata <- makejagsDRnetmeta(data=mydata,cov=F,class=F,ref='placebo')
# knots=c(20,30,50)
# dosresnetmeta.jagsdata$nn <- dosresnetmeta.jagsdata$n
# dosresnetmeta.jagsdata$rr <- dosresnetmeta.jagsdata$r
# dosresnetmeta.jagsdata$new.dose <- seq(1,80,1)
# dosresnetmeta.jagsdata$f.new.dose <- rcspline.eval(dosresnetmeta.jagsdata$new.dose,knots,inclx = T)[,2]
# dosresnetmeta.jagsdata$nd.new <- length(dosresnetmeta.jagsdata$new.dose)
#** Run JAGS
# 1. Run DRNMA - cubic splines
#dosresnetmeta.jagsdata$ref <- 15
# dosresnetmeta.jagsdata$trt1 <- unique(dosresnetmeta.jagsdata$t[,1])
# dosresnetmeta.jagsdata$trt2 <- unique(c(dosresnetmeta.jagsdata$t[,-1]))[!is.na(unique(c(dosresnetmeta.jagsdata$t[,-1])))]
dosresnetmeta.jagsdata$nd.new
dosresnetmeta.runjagsRCS <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('d','d2','dd','dd2','tau','Z','p.drug','cumeffectiveness'),model.file = modelDRnetmetaSpline,
n.chains=3,n.iter = 1000,n.burnin = 200,DIC=F,n.thin = 1)
dosresnetmeta.runjagsRCSInc <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('dd','dd2','tau','Z','p.drug','cumeffectiveness'),model.file = modelDRnetmetaSplineInc,
n.chains=3,n.iter = 1000,n.burnin = 200,DIC=F,n.thin = 1)
dosresnetmeta.runjagsRCS$BUGSoutput$summary
names.inc1 <- paste0('dd[',rep(dosresnetmeta.jagsdata$trt1,each=length(dosresnetmeta.jagsdata$trt2)),',',dosresnetmeta.jagsdata$trt2,']')
names.inc2 <- paste0('dd2[',rep(dosresnetmeta.jagsdata$trt1,each=length(dosresnetmeta.jagsdata$trt2)),',',dosresnetmeta.jagsdata$trt2,']')
est.con1 <- dosresnetmeta.runjagsRCS$BUGSoutput$summary[rownames(dosresnetmeta.runjagsRCS$BUGSoutput$summary)%in%names.inc1,]
est.con2 <- dosresnetmeta.runjagsRCS$BUGSoutput$summary[rownames(dosresnetmeta.runjagsRCS$BUGSoutput$summary)%in%names.inc2,]
est.incon1 <- dosresnetmeta.runjagsRCS$BUGSoutput$summary[names.inc1,]
est.incon2 <- dosresnetmeta.runjagsRCS$BUGSoutput$summary[names.inc2,]
plot(abs(est.con1[,'mean']),abs(est.incon1[,'mean']),ylim = c(0,0.03),xlim = c(0,0.03))
abline(0,1)
myd <- data.frame(y=est.con2[,'mean'],x=est.incon2[,'mean'])
ggplot(myd,aes(y=y,x=x))+
geom_point()+
geom_abline(slope=1,intercept=0)+
theme_bw()
# 2. Run DRNMR - cubic splines - ROB
mydata$covariate <- as.numeric(as.factor(antidepdata$Overall_study_RoB))-2
dosresnetmeta.jagsdata <- makejagsDRnetmeta(data=mydata,cov=T)
dosresnetmeta.jagsdata$ref <- 5
dosresnetmetaregROB.runjagsRCS <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('d','d2','bb','tua'),model.file = modelDRnetmetaregSpline,
n.chains=3,n.iter = 1000,n.burnin = 200,DIC=F,n.thin = 1)
dosresnetmetaregROB.runjagsRCS
# 3. Run DRNMR - cubic splines - VAR: small-study effect
mydata$covariate <- antidepdata$selogOR^2
dosresnetmeta.jagsdata <- makejagsDRnetmeta(data=mydata,cov=T)
dosresnetmeta.jagsdata$ref <- 5
dosresnetmetaregVAR.runjagsRCS <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('d','d2','bb','tau'),model.file = modelDRnetmetaregSpline,
n.chains=3,n.iter = 1000,n.burnin = 200,DIC=F,n.thin = 1)
dosresnetmetaregVAR.runjagsRCS
# 4. Run DRNMR - cubic splines - year of study **** Does not run, study_year has NA's
mydata$covariate <- antidepdata$Study_year-min(antidepdata$Study_year,na.rm=T)
dosresnetmeta.jagsdata <- makejagsDRnetmeta(data=mydata[!is.na(mydata$covariate),],cov=T)
dosresnetmeta.jagsdata$ref <- 5
dosresnetmetaregYEAR.runjagsRCS <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('d','d2','bb','tau'),model.file = modelDRnetmetaregSpline,
n.chains=3,n.iter = 1000,n.burnin = 200,DIC=F,n.thin = 1)
dosresnetmetaregYEAR.runjagsRCS
# 5. Run DRNMA - cubic splines + class effect
mydata$class <- class
dosresnetmeta.jagsdata <- makejagsDRnetmeta(data=mydata,cov=F,includeClass=T)
dosresnetmeta.jagsdata$ref <- 5
dosresnetmetaClass.runjagsRCS <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('dd','dd2','sd'),model.file = modelDRnetmetaSplineClass,
n.chains=3,n.iter = 1000,n.burnin = 200,DIC=F,n.thin = 1)
dosresnetmetaClass.runjagsRCS
# 6. Run DRNMR - cubic splines + class effect + covariate - ROB
mydata$covariate <- as.numeric(as.factor(antidepdata$Overall_study_RoB))-2
dosresnetmeta.jagsdata <- makejagsDRnetmeta(data=mydata,cov=T,includeClass=T)
dosresnetmeta.jagsdata$nc <- 3
dosresnetmeta.jagsdata$ref <- 5
dosresnetmetaregClass.runjagsRCS <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('dd','dd2','sd'),model.file = modelDRnetmetaregSplineClass,
n.chains=3,n.iter = 1000,n.burnin = 200,DIC=F,n.thin = 1)
dosresnetmetaregClass.runjagsRCS
#
as.numeric(as.factor(antidepdata$Drug))
#
sort(levels(as.factor(antidepdata$Drug)))
#
sort(c('c','b','a'))
install_github("audrey-b/BUGSnet",force=TRUE)
library(BUGSnet)
data("thrombolytic")
head(thrombolytic)
myd <- data.prep(antidep,varname.t = 'drug',varname.s = 'studyid')
net.plot(myd,node.scale = 0.2,edge.scale = 0.25)
#
# Fractional polynomials
library(mfp)
library(dosresmeta)
str(antidepdata)
f <- mfp(logRR~fp(hayasaka_ddd),data = antidepdata,family = gaussian)
summary(f)
plot(f)
antidepdata$hayasaka_ddd2 <- ifelse(antidepdata$hayasaka_ddd ==0, 0,antidepdata$hayasaka_ddd^-3)
ff <- dosresmeta(formula=logRR~I(hayasaka_ddd)+hayasaka_ddd2, proc="1stage",id=Study_No, type=type,cases=Responders,n=No_randomised,se=selogRR,data=antidepdata,method = 'reml')
ff2 <- dosresmeta(formula=logRR~rcs(hayasaka_ddd,knots), proc="1stage",id=Study_No, type=type,cases=Responders,n=No_randomised,se=selogRR,data=antidepdata,method = 'reml')
#
knots=c(20,30,50)
newdose <- 1:80
dd <- 0.0058*newdose+10.4571*newdose^-2
dd1 <- 65.8204*newdose^-2
dd3 <- 0.0001*newdose^2+46.9503*newdose^-2
dd2 <- 0.0091*newdose+-0.0156*rcspline.eval(newdose,knots,inclx = T)[,2]
with(antidepdata,plot(hayasaka_ddd,logRR))
with(antidepdata,lines(dd))
with(antidepdata,lines(dd1,col=3))
with(antidepdata,lines(dd2,col=2))
with(antidepdata,lines(dd3,col=4))
#
# predict the absolute effect + ranking: SUCRA
# dosresnetmeta.runjagsRCS <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('d','d2','Z','sd','p.drug','SUCRA'),model.file = modelDRnetmetaSpline,
# n.chains=3,n.iter = 1000,n.burnin = 200,DIC=F,n.thin = 1)
# dosresnetmeta.runjagsRCS$BUGSoutput$mean$SUCRA
#
# p.placebo<- exp(dosresnetmeta.runjagsRCS$BUGSoutput$mean$Z)/(1+exp(dosresnetmeta.runjagsRCS$BUGSoutput$mean$Z))
# p.drug <- dosresnetmeta.runjagsRCS$BUGSoutput$mean$p.drug
# l.ci <- dosresnetmeta.runjagsRCS$BUGSoutput$summary[paste0('p.drug[',1:80,',',6,']'),'2.5%']
# u.ci <- dosresnetmeta.runjagsRCS$BUGSoutput$summary[paste0('p.drug[',1:80,',',6,']'),'97.5%']
# lp.ci <- exp(dosresnetmeta.runjagsRCS$BUGSoutput$summary['Z','2.5%'])/(1+exp(dosresnetmeta.runjagsRCS$BUGSoutput$summary['Z','2.5%']))
# up.ci <- exp(dosresnetmeta.runjagsRCS$BUGSoutput$summary['Z','97.5%'])/(1+exp(dosresnetmeta.runjagsRCS$BUGSoutput$summary['Z','97.5%']))
#
# library(ggplot2)
# df1 <- data.frame(new.dose=dosresnetmeta.jagsdata$new.dose,y1 =p.drug[,5],y2=l.ci,
# y3=u.ci,yp1=p.placebo,yp2=lp.ci,yp3=up.ci)
#
# ggplot(data = df1,aes(x=new.dose)) +
# geom_line(aes(y=y1,color='treatment'))+
# geom_line(aes(y=yp1,color='placebo'))+
# scale_color_manual(values=c('steelblue','coral4'))+
# geom_smooth(aes(x=new.dose, y=y1, ymin=y2,
# ymax=y3),color='coral4',fill='coral1',
# data=df1, stat="identity")+
# geom_smooth(aes(x=new.dose, y=yp1, ymin=yp2,
# ymax=yp3),color='steelblue',fill='steelblue2',
# data=df1, stat="identity")+
# xlab('')+
# ylab('')+
# ylim(0,0.9)+
# theme(panel.background = element_rect(fill = 'snow1',colour = 'white'),
# legend.position = c(0.72,0.95), legend.key.size = unit(3, "cm"), legend.text = element_text(size=12),
# legend.key.height = unit(0.5,"cm"),legend.title = element_blank(), legend.background = element_blank(),
# axis.text.x = element_text(face='bold',size=14),
# axis.text.y = element_text(face='bold',size=14))
# Timesheet_Template_2019 <- read_excel("~/Google Drive/Administrative/Timesheet_Template_2019.xlsx",
# sheet = "JAN", col_names = FALSE)
# Timesheet_Template_2019 <- as.data.frame(Timesheet_Template_2019)
# Timesheet_Template_2019[c('Work Package 1'),2:(ncol(Timesheet_Template_2019)-2)] <- sample(seq(7,9,by =0.5),ncol(Timesheet_Template_2019)-3,replace = T)
# for (i in dosresnetmeta.jagsdata$trt2) {
# print(i)
# }
#
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.